#include "../include/Multigrid_Solver.h"
#include "../include/Multigrid_Solver_1D.h"
#include "../include/Multigrid_Solver_2D.h"
#include "../include/Multigrid_Factory.h"
#include "../include/Functions_test.cpp"


int main()
{
    Factory1DMGSolver Solver1D;
    Factory2DMGSolver Solver2D;
    
    FactoryMGSolverBase *_p_fac_base = &Solver1D;
    Multigrid_Solver * _arena1D = _p_fac_base ->CreateMGSolver();
    _p_fac_base = &Solver2D;
    Multigrid_Solver * _arena2D = _p_fac_base->CreateMGSolver();

double _a[1]={0};
int n=128;
    VectorXd u_,f_,r_,r_origin,acu;
    acu.setZero(n+1);
for(int i=0;i<=n;i++)
{   
    _a[0]=1.0*i/n;
    acu(i)=Dim1_type1_f(_a);
};

int _Rtype;
int _Itype;
std::cout<<"Give me n and Rtype and Itype,:"<<std::endl;
std::cin>>n;
std::cin>>_Rtype;
std::cin>>_Itype;

_arena2D->Initialize(n,Dim2_type1_df2rd,Dim2_type1_f);
_arena2D->setRestircType(_Rtype);
_arena2D->setIntpType(_Itype);
_arena2D->setSmoothType(2);

double acu_norm=acu.lpNorm<Eigen::Infinity>();
double h;
double r_origin_norm;
int cnt;
u_=_arena2D->get_u_origin();
f_=_arena2D->get_f_origin();
h = _arena2D->get_h();
r_origin=_arena2D->ResidualSolver(u_,f_,h);
r_origin_norm=r_origin.lpNorm<Eigen::Infinity>();

    u_=_arena2D->FMGcycle(f_,h);
        r_=_arena2D->ResidualSolver(u_,f_,h);   
std::cout<<"n= "<<n<<std::endl;
std::cout<<"Residual norm_inf="<<r_.lpNorm<Eigen::Infinity>()<<std::endl;

std::cout<<"relative accu= "<<r_.lpNorm<Eigen::Infinity>() /r_origin_norm<<std::endl;

// u_=_arena2D->get_u_origin();
// for (cnt=1;cnt<=100;cnt++)
// {
//     u_=_arena2D->Vcycle(u_,f_,h);
//     r_=_arena2D->ResidualSolver(u_,f_,h);   
//     if (r_.lpNorm<Eigen::Infinity>() <= r_origin_norm*1e-8)
//     break;
// };

// std::cout<<"n= "<<n<<std::endl;
//  std::cout<<"cnt= "<<cnt<<std::endl;
// //  std::cout<<wucha<<std::endl;
// std::cout<<"Residual norm_inf="<<r_.lpNorm<Eigen::Infinity>()<<std::endl;

// std::cout<<"relative accu= "<<r_.lpNorm<Eigen::Infinity>() /r_origin_norm<<std::endl;

_arena1D->Initialize(n,Dim1_type1_df2rd,Dim1_type1_f);
_arena1D->setRestircType(_Rtype);
_arena1D->setIntpType(_Itype);
_arena1D->setSmoothType(2);


u_=_arena1D->get_u_origin();
f_=_arena1D->get_f_origin();
h = _arena1D->get_h();
r_origin=_arena1D->ResidualSolver(u_,f_,h);
r_origin_norm=r_origin.lpNorm<Eigen::Infinity>();

for (cnt=1;cnt<=100;cnt++)
{
    u_=_arena1D->Vcycle(u_,f_,h);
    r_=_arena1D->ResidualSolver(u_,f_,h);   
    if (r_.lpNorm<Eigen::Infinity>() <= r_origin_norm*1e-8)
    break;
};

std::cout<<"n= "<<n<<std::endl;
 std::cout<<"cnt= "<<cnt<<std::endl;
//  std::cout<<wucha<<std::endl;
std::cout<<"Residual norm_inf="<<r_.lpNorm<Eigen::Infinity>()<<std::endl;

std::cout<<"relative accu= "<<r_.lpNorm<Eigen::Infinity>() /r_origin_norm<<std::endl;
    
return 0;
};

